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Method and System for Numerically Simulating Foam-like Material in 

Finite Element Analysis 

CROSS REFERENCE TO RELATED APPLICATION 

[0001] This application is a continuation-in-part of a U.S. patent application 
Ser. No. 10/669,047 filed Sep. 23, 2003, the entire contents of which are 
incorporated herein by reference. 

BACKGROUND OF THE INVENTION 

1. Field of the Invention 

[0002] The present invention generally relates to method, system and 
software product used in finite element analysis, more particularly to numerical 
simulation of the structural behaviors of highly compressible material such as 
foam under loading conditions. 

2. Description of the Related Art 

[0003] Finite element analysis (FEA) is a computerized method widely used in 
industry to model and solve engineering problems relating to complex systems. 
FEA derives its name from the manner in which the geometry of the object 
under consideration is specified. With the advent of the modern digital 
computer, FEA has been implemented as FEA software. Basically, the FEA 
software is provided with a model of the geometric description and the 
associated material properties at each point within the model. In this model, the 
geometry of the system under analysis is represented by solids, shells and 
beams of various sizes, which are called elements. The vertices of the elements 
are referred to as nodes. The model is comprised of a finite number of 
elements, which are assigned a material name to associate with material 
properties. The model thus represents the physical space occupied by the 
object under analysis along with its immediate surroundings. The FEA software 
then refers to a table in which the properties (e.g., stress-strain constitutive 
equation, Young's modulus, Poisson's ratio, thermo-conductivity) of each 
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material type are tabulated. Additionally, the conditions at the boundary of the 
object (i.e., loadings, physical constraints, etc.) are specified. In this fashion a 
model of the object and its environment is created. 

[0004] FEA is becoming increasingly popular with automobile manufacturers 
for optimizing both the aerodynamic performance and structural integrity of 
vehicles. Similarly, aircraft manufacturers rely upon FEA to predict airplane 
performance long before the first prototype is built. Rational design of 
semiconductor electronic devices is possible with Finite Element Analysis of the 
electrodynamics, diffusion, and thermodynamics involved in this situation. FEA 
is utilized to characterize ocean currents and distribution of contaminants. FEA 
is being applied increasingly to analysis of the production and performance of 
such consumer goods as ovens, blenders, lighting facilities and many plastic 
products. In fact, FEA has been employed in as many diverse fields as can be 
brought to mind, including plastics mold design, modeling of nuclear reactors, 
analysis of the spot welding process, microwave antenna design, simulating of 
car crash and biomedical applications such as the design of prosthetic limbs. In 
short, FEA is utilized to expedite design, maximize productivity and efficiency, 
and optimize product performance in virtually every stratum of light and heavy 
industry. This often occurs long before the first prototype is ever developed. 
[0005] The finite element analysis method is described in detail by Thomas J. 
R. Hughes in 'The Finite Element Method" (1987), published by Prentice-Hall, 
Inc., New Jersey, which is incorporated herein by reference in its entirety. 
Generally, FEA begins by generating a finite element model of a system. In this 
model, a subject structure is reduced into a number of node points which are 
connected together to form finite elements. Governing equations of motion are 
written in a discrete form, where the motions of each node point are the 
unknown part of the solution. A simulated load or other influence is applied to 
the system and the resulting effect is analyzed using well known mathematical 
methods. 

[0006] FEA software can be classified into two general types, implicit FEA 
software and explicit FEA software. Implicit FEA software uses an implicit 
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equation solver to solve a system of coupled linear equations. Such software is 
generally used to simulate static or quasi-static problems. Explicit FEA 
software does not solve coupled equations but explicitly solves for each 
unknown assuming them uncoupled. Explicit FEA software usually uses central 
difference time integration which requires very small solution cycles or time 
steps for the method to be stable and accurate. Explicit FEA software is 
generally used to simulate short duration events where dynamics are important 
such as impact type events. 

[0007] One of the most challenging FEA tasks is to simulate an impact event 
such as car crash. The highly non-linear behavior of the structural materials 
must be numerically simulated accurately, realistically and efficiently. A number 
of materials such as steel, aluminum, foam and rubber that used in an 
automobile must be included in such FEA software. Many components used in 
an automobile are made of materials such as foam and rubber, the simulation 
of the structural responses of these materials becomes very important for the 
overall accuracy of an analysis. 

[0008] Traditionally the structural responses of highly compressible material 
such as foam have been numerically simulated using the Ogden strain energy 
function W. A brief summary of the Ogden strain energy function and 
corresponding engineering/nominal stress ao and Cauthy/true stress a functions 
are listed in FIG. 1 . More details of the Ogden energy function is described by 
R.W. Ogden in Chapter 7 of the book titled: "Non-linear Elastic Deformations" 
(1984), published by Ellis Horwood Limited, United Kingdom, which is 
incorporated herein by reference in its entirety. A number of commercially 
available FEA software includes these approaches to simulate foam-like 
material. For example, LS-DYNA, a general purpose three-dimensional non- 
linear large-deformation FEA software from Livermore Software Technology 
Corporation, is capable of simulating foam-like material using the Ogden energy 
function. 

[0009] Today, there are a number of practical problems associated with the 
simulation of foam-like material in FEA. To implement the Ogden energy 
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function properly in the FEA software requires engineers to spend a 
tremendous amount of time to prepare experimental data and then convert 
them into a set of coefficients to fit a polynomial Ogden function for FEA 
software. Due to highly non-linear characteristics of this polynomial function, 
the inexactly fitted function has often resulted. This leads to a lengthy iterative 
trial-and-error process of modifying the input coefficients to match the behavior 
of foam-like material in the real world. It is therefore desirable to have a new 
method to numerically simulate foam-like material more efficiently and 
effectively. 

SUMMARY OF THE INVENTION 

[0010] This section is for the purpose of summarizing some aspects of the 
present invention and to briefly introduce some preferred embodiments. 
Simplifications or omissions may be made to avoid obscuring the purpose of the 
section. Such simplifications or omissions are not intended to limit the scope of 
the present invention. 

[0011] The present invention discloses a system, method and software 
product to numerically simulate highly compressible material such as foam with 
a pragmatic approach. Under the assumption of uniaxial loading and isotropic 
material, a method of calculating stress function f{X) is developed as shown is 
FIGS. 2 and 3. The method eliminates the requirement of fitting a polynomial 
function; instead f(A,) is calculated using only a set of the stress-strain curve 
from a uniaxial loading test. Based upon the data of this stress-strain curve, the 
stress function f(k) is pre-computed for a range of equally spaced stretch ratios 
X. These results are stored in a computer's RAM as a lookup table and retrieved 
later during the solution phase of FEA. The present invention does not require 
time-consuming trial-and-error process of fitting polynomial coefficients. This is 
a huge advantage over the exiting method. The stresses corresponding to a 
particular strain or stretch ratio can be interpolated rapidly from the pre- 
computed table directly during the lengthy solution phase. 
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[0012] According to one aspect of the present invention, during the solution 
phase of FEA, the stresses at each integration point of the finite elements is 
computed for each time step. In order to compute the stresses, the 
eigensolution of a 3x3 stretch tensor is computed first. The resulting 
eigenvalues are the three main stretch ratios for the element. The 
corresponding engineering/nominal stress is then looked up or interpolated. 
These nominal stresses are converted into the true stresses and transformed to 
the global coordinate system. If the rate dependent viscoelastic effect is 
needed, the nominal stresses are looked up from a set of curves instead of a 
single curve. 

[0013] Other objects, features, and advantages of the present invention will 
become apparent upon examining the following detailed description of an 
embodiment thereof, taken in conjunction with the attached drawings. 



BRIEF DESCRIPTION OF THE DRAWINGS 

[0014] These and other features, aspects, and advantages of the present 
invention will be better understood with regard to the following description, 
appended claims, and accompanying drawings as follows: 

[0015] FIG. 1 lists the Ogden strain energy function and associated formula for 
nominal and true stresses for compressible material such as foam. 

[0016] FIG. 2 and 3 show the detailed derivation of the stress function f(A,) for 
the present invention based on compressible material and uniaxial loading. The 
stress function f(?i) is only dependent upon stress-strain curve obtain via a 
uniaxial loading test. 

[0017] FIG. 4 shows an exemplary engineering stress-strain curve obtained 
from a uniaxial loading test. 

[0018] FIG. 5 shows a process to calculate the stress function f(k) based on 
the present invention. 
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[0019] FIG. 6 is an exemplary flow chart of an implementation of the present 
invention in FEA software. 

[0020] FIG. 7 depicts a block diagram of an exemplary computer, in which the 
present invention may be implemented. 
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DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS 

[0021] In the following description, numerous specific details are set forth in 
order to provide a thorough understanding of the present invention. However, it 
will become obvious to those skilled in the art that the present invention may be 
practiced without these specific details. The descriptions and representations 
herein are the common means used by those experienced or skilled in the art to 
most effectively convey the substance of their work to others skilled in the art. 
In other instances, well-known methods, procedures, components, and circuitry 
have not been described in detail to avoid unnecessarily obscuring aspects of 
the present invention. 

[0022] Reference herein to "one embodiment" or "an embodiment" means that 
a particular feature, structure, or characteristic described in connection with the 
embodiment can be included in at least one embodiment of the invention. The 
appearances of the phrase "in one embodiment" in various places in the 
specification are not necessarily all referring to the same embodiment, nor are 
separate or alternative embodiments mutually exclusive of other embodiments- 
Further, the order of blocks in process flowcharts or diagrams representing one 
or more embodiments of the invention do not inherently indicate any particular 
order nor imply any limitations in the invention. 
[0023] To facilitate the description of the present invention, it deems 
necessary to provide definitions for some terms that will be used throughout the 
disclosure herein. It should be noted that the definitions following are to 
facilitate the understanding and describe the present invention according to an 
embodiment. The definitions may appear to include some limitations with 
respect to the embodiment, the actual meaning of the terms has applicability 
well beyond such embodiment, which can be appreciated by those skilled in the 
art: 

[0024] FEA stands for Finite Element Analysis. 

[0025] Implicit FEA refers to K x = F, where K is the global stiffness matrix, x 
is the unknown displacement array and F is the global force array. 
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[0026] Explicit FEA refers to M a = F, where M is the mass array, a is the 
acceleration array and F is the global force array. 
[0027] Time domain analysis refers to a FEA that simulates a physical 
phenomenon in time domain. 

[0028] Dynamic analysis refers to a FEA that accounts for the mass and 
inertia effects of a structure. In general, there are two classes of dynamic 
analysis: time domain and frequency domain. 

[0029] Solution cycle and cycle are used interchangeably to refer to the 
numbered solution states starting with cycle 0 at time 0. 
[0030] The time step refers to an interval of time between subsequent cycles. 
[0031] Loading Condition is defined as the external load acting on a structure. 
[0032] Uniaxial Load is either tension or compression in one direction. 
[0033] Longitudinal Direction is the direction, in which the uniaxial force 
applies to a material. 

[0034] Lateral Direction is the orthogonal direction to the longitudinal direction. 

[0035] Tension is a force to pull or stretch a material. 

[0036] Compression is a force to squeeze or press a material. 

[0037] Stress is defined as the intensity of force or force per unit area of the 

material. 

[0038] Stretch Ratio X is defined as the ratio between the deformed length 
and the original length of a material. 

[0039] Strain e is a non-dimensional quantity representing the elongation or 
shrinkage per unit length of a material under loading condition. The Strain and 
Stretch Ratio is related with the following equation: s = . 
[0040] Strain Rate is defined as the amount of strain that develops in one unit 
of time. 

[0041] Poisson's Ratio v is defined as a function of the strain in lateral 
direction to the strain in the longitudinal direction. Strain Energy is potential 
energy generated from the force acting on a material. In the elastic case, the 
strain energy is defined as W = 0.5*e*a. 
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[0042] Deformation Gradient Matrix/Tensor is defined as the local nature of 
the deformation of a solid body. 

[0043] Principal Stretch Ratio is defined as the stretch ratio in the principal 
loading direction. 

[0044] Eigensolution is a mathematical method to determine the natural 

vibration mode shapes and amplitudes of a structure. 

[0045] Eigenvalues are the results obtained from the eigensolution. 

[0046] Elastic material is a material recovers its original shape after 

unloading. 

[0047] Hyperelastic material is a material has one-to-one relationship between 
stress and strain for a large deformation. The stress value can be determined 
by an energy function such as Stress function f(?i). 

[0048] Quasistatic loading is defined as the static force dominates comparing 
to inertia and damping forces. 

[0049] Rate dependent under dynamic loading means that under dynamic 
loading the material is no longer quasistatic; for the same values of strain, 
different values of stress may be obtained depending upon the loading velocity 
or strain rate. Rate dependent under dynamic loading and viscous material are 
interchangeably used. 

[0050] Nominal stress or Engineering stress, interchangeably used herein, is 

defined as the force divided by undeformed cross section area. 

[0051] True stress or Cauchy stress, interchangeably used herein, is defined 

as the force divided by current or deformed cross section area. 

[0052] Local coordinate system or element coordinate system, 

interchangeably used herein, is a coordinate system used to define element 

structural responses such as true stress and nominal stress. 

[0053] Global coordinate system is a coordinate system used to define the 

structure in FEA. 

[0054] Embodiments of the present invention are discussed herein with 
reference to FIGS. 1-7. However, those skilled in the art will readily appreciate 
that the detailed description given herein with respect to these figures is for 



10 



explanatory purposes as the invention extends beyond these limited 
embodiments. 

[0055] Referring now to the drawings, in which like numerals refer to like parts 
throughout several views. FIG. 1 shows the structural response functions in the 
original Ogden highly compressible material model and derivation of nominal 
stress formula as function of longitudinal stretch ratio h based on the 
assumption of uniaxial loading and isotropic material. Ogden strain energy 
function W for compressible material is listed in 110. In these functions, h, ^2 
and fa represent the stretch ratios in three principal directions of a compressible 
material. The variable J is the volumetric ratio or the deformed volume divided 
by the undeformed volume. The material coefficient n is related to Poisson's 
ratio of the material. Coefficients \i and a are the material constants. The 
classic strain energy function W in 110 is the foundation to numerically simulate 
highly compressible material using a digital computer. The formula for true 
stress o is listed in 120 and nominal stress ao is listed in 130. Under the 
uniaxial tension or compressible loading, the coefficient n is determined in the 
formulas listed in 140. Under the isotropic material assumption, the two lateral 
stretch ratios X2 and X3 are equal to each other. The relationship between 
stretch ratio X.1 in the uniaxial loading direction and lateral stretch ratio X3 is also 
listed in 140. Substitute X 2 and A. 3 with in the nominal stress formula in 130, 
the nominal stress formula can be written as a function of the longitudinal 
stretch ratio X.1 as shown in 150. It is evident that element stresses can easily be 
calculated, once the coefficients n and a have been determined; however, the 
determination of these coefficients for the polynomial function is very difficult in 
practice. 

[0056] FIGS. 2 and 3 in combination show the derivation of stress function 
f(A,). The original stress function f(A.) in polynomial form is defined in 210. 
Substitute into the nominal stress a 0 formula in 150, the resulting formula is 
shown in 220. In 230, Poisson's ratio v is defined as function of material 
constant n. After substituting v in formula 220, the new formula is listed as 
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formula 240. Formula 250 is a result of replacing X with X' v in formula 240. 
The next formula 260 is similarly obtained by substituting a ( ~ u)2 into X in formula 
240. The process repeats for the rest of higher order formulas using the same 
substitution method. 

[0057] By rearranging formulas 240, 250, 260 and the rest of higher order 
formulas, the resulting equation for stress function f(X) is rewritten as listed in 
320. The equation 320 is rewritten as a form equation as shown in 330. 
Wherein the value of stress function i{X) for the stretch ratio of interest X equals 
to the summation of a sequence of x { ~ v)} ^(^J , where j is an integer related to 

the j-th term of the sequence, v is Poisson's ratio of the compressible material, 
and is the stress value at stretch ratio x { ~ v>J . Using the relationship 

between strain e and stretch ratio X in 340, the equation 330 transforms to 
equation 350, which can be evaluated solely with the stress values 
corresponding to the strains from the strain-stress curve obtained via a uniaxial 
tension/compression test of the material of interest. Finally, the nominal stress 
go and true stresses a are rewritten and dependent only on the stress function 
f(X) listed in formulas 350 and 360, respectively. 

[0058] Referring now to FIG. 4, an engineering or nominal stress-strain curve 
420 is illustrated. The curve represents a set of experimental data for a 
compressible material under uniaxial loading condition; both tension and 
compression tests results are included. According to one embodiment of the 
present invention, the engineer needs only to prepare the experimental data 
such as this exemplary curve 420 as input to simulate compressible material in 
FEA software. In order to capture the stretch ratio of interest, an exemplary 
input requires a range of strain values shown in 430. If the viscous effect of 
compressible material is under consideration, the exemplary curve 420 may be 
replaced by a family of stress-strain curves resulting from dynamic experiments. 
[0059] With reference now to FIG. 5, it shows a flow chart 500 for computing 
the infinite series in equation 330. Because the order of magnitude decreases 
drastically from one term to the next, the sequence converges very rapidly. 
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According to one embodiment of the present invention, the flow chart 500 of 
evaluation of function i(X) is summarized in FIG. 5. At 510, f(A.) is assigned a 
value equals to X multiplied by oo(A,-1) for a given X. The stress oo value at 
strain eoor A.-1 is from engineering test strain-stress data (e.g., the curve 
defined in the user input phase of FEA software). At 520, X is stored into a 
variable W A new variable X n&N is stored as W v at 530, where v is Poisson's 
ratio of the compressible material. At 540, a comparing test is performed for the 
absolute value of (>,new-1 ) being less than or equal to a threshold value £. In 
one embodiment, £ is set to 0.01 . If the test succeeds, the computation of f(A.) 
has finished, the rest of the terms in the infinite series is too small to affect the 
final result of the computation. If the test fails, another stress value at strain 
Xnew-1 is multiplied by X n e W and accumulated into function 1(X) at 550. At 560, 
the value of X nevi is stored into Xow- The process goes back to 530 until the 
computation finishes. 

[0060] An exemplary FEA software implementation 600 of the present 
invention is shown in FIG. 6. At the input phase 610 of FEA software, at least 
one set of engineering or nominal stress-strain data for rubber material under 
uniaxial loading test is read in as a group of corresponding stress-strain pairs. If 
the loading condition is quasistatic, only one set of stress-strain curve is 
required. If the rate dependent dynamic loading condition is under 
consideration, then a family of stress-strain curves representing different strain 
rate are required. Based on the input stress-strain constitutive curve, a plurality 
of f{X) is then generated for a range of X at regular intervals at the initialization 
phase 620 using the process listed in flow chart 500 as shown in FIG. 5. These 
pre-computed f(A,) are stored in computer's RAM as a stress function table for 
later use in the solution phase. Next in 630, the FEA software calculates 
structural responses at each integration point of every element for each time 
step. According to one embodiment, for the element with compressible material 
such as foam, a subroutine handling foam material is called. Within the 
subroutine, the deformation gradient matrix is gathered first at 640. For a solid 
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element, either the left or the right 3x3 stretch tensors may represent the 
deformation gradient matrix. The deformation gradient matrix may be 
corresponded to the stretch ratio X or the stretch ratio squared X 2 depending 
upon the solution scheme employed. At 650, an eigensolution is solved for the 
deformation gradient tensor. The principal stretch ratio can be calculated from 
the resulting eigenvalues at 660. With these stretch ratios, at 670, the nominal 
stresses are then calculated using the corresponding 1(X) values stored in the 
initialization phase. Because all f(X) values of interest are pre-computed and 
stored as tabulated form, it is very computational efficient to perform a table 
lookup and interpolation. In the case of rate dependent dynamic loading, strain 
rate at the element integration point is used as basis to select two closest 
curves to determine the corresponding rate dependent stress. The true 
stresses are calculated and transformed into the global coordinate system at 
690. 

[0061] With reference now to FIG. 7, a block diagram illustrates a computing 
device 700 in which the present invention may be implemented, and in which 
code or instructions implementing the processes of the present invention may 
be located. The exemplary computer system in FIG. 7 is discussed only for 
descriptive purposes. It should not be considered a limitation of the invention. 
Although the following descriptions related to a particular computer system, the 
concepts apply equally to other computer systems that are dissimilar to that 
shown in FIG. 7. 

[0062] Computer system 700 includes a processor 710 and main random 
access memory (RAM) 720 connecting to a local bus 705 through a bridge 715. 
Additional connections to local bus 705 may be made through direct component 
interconnection or through add-in boards. In the depicted example, network 
adapter 725, small computer system interface (SCSI) adapter 730, and 
expansion bus interface 735 are directly connected to local bus 705. In contrast, 
audio adapter 740, graphics adapter 745, and video adapter 750 are connected 
to local bus 705 by add-in boards inserted into expansion slots. Expansion bus 
interface 735 provides a connection for a keyboard and mouse adapter 755, 
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modem 760, and additional memory 765. SCSI adapter 730 provides a 
connection for hard disk drive 770, tape drive 775, and CD-ROM drive 780. 
[0063] In order to communicate with other computer systems via a network, 
the computer system 700 connects to the network via network adapter 725. 
The network, Internet or intranet, connects multiple network devices utilizing 
general purpose communication lines. 

[0064] Those of ordinary skill in the art will appreciate that the hardware 
shown in FIG. 7 may vary depending on the implementation. Other internal 
hardware or peripheral devices, such as flash ROM (or equivalent nonvolatile 
memory) or optical disk drives and the like, may be used in addition to or in lieu 
of the hardware depicted in FIG. 7. Also, the processes of the present invention 
may be applied to a multiprocessor computer system. In general, Computer 
system 700 is controlled and coordinated by operating system (OS) software, 
which performs tasks such as process scheduling, memory management, 
networking and I/O services. Exemplary OS includes Linux™, Microsoft 
Windows™. 

[0065] Although an exemplary embodiment of invention has been disclosed, it 
will be apparent to those skilled in the art that various changes and 
modifications may be made to achieve the advantage of the invention. It will be 
obvious to those skilled in the art that some components may be substituted 
with another component providing same function. The appended claims cover 
the present invention. 
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